# Load PRL theme
source("Code/ggplot2_theme.R")

# Load bootstrap results
boot_results <- readRDS("Data/panel/bootstrapped_within_ind_SDs.rds")

# Load global means
global_means <- readRDS("Data/panel/global_means.rds")

####################
### PREPARE DATA ###
####################

# Rename variables
global_means <- global_means %>%
    rename(Variable = skim_variable) %>%
    mutate(
        Variable = case_when(
            Variable == "in_party_therm" ~ "Warmth\nToward In-Party",
            Variable == "norm_censorship" ~ "Belief that Government\nShould Censor Out-Party Media",
            Variable == "norm_censorship_perception" ~ "Perceived Out-Party Belief that\nGovernment Should Censor In-Party Media",
            Variable == "norm_companies" ~ "Belief that Officials\nShould Punish Out-Party Companies",
            Variable == "norm_companies_perception" ~ "Perceived Out-Party Belief that\nOfficials Should Punish In-Party Companies",
            Variable == "norm_executive" ~ "Belief that Presidents Should\nCircumvent Out-Party Congresses",
            Variable == "norm_executive_perception" ~ "Perceived Out-Party Belief that\nPresidents Should Circumvent In-Party Congresses",
            Variable == "norm_judges" ~ "Belief that Official Should\nIgnore Out-Party Court Decisions",
            Variable == "norm_judges_perception" ~ "Perceived Out-Party Belief that\nOfficials Should Ignore In-Party Court Decisions",
            Variable == "norm_loyalty" ~ "Belief that In-Partisans Should be Loyal to\nIn-Party Candidates Who Question Elections",
            Variable == "norm_loyalty_perception" ~ "Perceived Out-Party Belief that Out-Partisans Should\nbe Loyal to Out-Party Candidates Who Question Elections",
            Variable == "norm_polling" ~ "Belief that In-Party Should\nReduce Polling Stations in Out-Party Areas",
            Variable == "norm_polling_perception" ~ "Perceived Out-Party Belief that Out-Party\nShould Reduce Polling Stations in In-Party Areas",
            Variable == "out_party_therm" ~ "Warmth\nToward Out-Party",
            Variable == "violence_protest" ~ "Support for In-Partisan Who Led\nan Illegal Protest Against the Out-Party",
            Variable == "violence_vandalism" ~ "Support for In-Partisan\nWho Vandalized Out-Party Signs",
            Variable == "violence_assault" ~ "Support for In-Partisan Who\nThrew Rocks at Out-Partisans",
            Variable == "violence_assault_perception" ~ "Perceived Out-Party Support for\nOut-Partisan Who Threw Rocks at In-Partisans",
            Variable == "violence_arson" ~ "Support for In-Partisan Who\nSet Fire to Out-Party Headquarters",
            Variable == "violence_assault_deadly_weapon" ~ "Support for In-Partisan\nWho Drove into Out-Party Crowd",
            Variable == "violence_murder" ~ "Support for In-Partisan Who\nStabbed Prominent Out-Partisan",
            Variable == "violence_murder_perception" ~ "Perceived Out-Party Support for\nIn-Partisan Who Stabbed Prominent Out-Partisan"
        )
    )

boot_results <- boot_results %>%
    pivot_longer(
        cols = -simulation,
        names_to = "Variable"
    ) %>%
    group_by(Variable) %>%
    reframe(
        estimate = value %>% mean(na.rm = TRUE),
        lower = value %>% quantile(.025),
        upper = value %>% quantile(.975)
    ) %>%
    mutate(
        Variable = case_when(
            Variable == "in_party_therm" ~ "Warmth\nToward In-Party",
            Variable == "norm_censorship" ~ "Belief that Government\nShould Censor Out-Party Media",
            Variable == "norm_censorship_perception" ~ "Perceived Out-Party Belief that\nGovernment Should Censor In-Party Media",
            Variable == "norm_companies" ~ "Belief that Officials\nShould Punish Out-Party Companies",
            Variable == "norm_companies_perception" ~ "Perceived Out-Party Belief that\nOfficials Should Punish In-Party Companies",
            Variable == "norm_executive" ~ "Belief that Presidents Should\nCircumvent Out-Party Congresses",
            Variable == "norm_executive_perception" ~ "Perceived Out-Party Belief that\nPresidents Should Circumvent In-Party Congresses",
            Variable == "norm_judges" ~ "Belief that Official Should\nIgnore Out-Party Court Decisions",
            Variable == "norm_judges_perception" ~ "Perceived Out-Party Belief that\nOfficials Should Ignore In-Party Court Decisions",
            Variable == "norm_loyalty" ~ "Belief that In-Partisans Should be Loyal to\nIn-Party Candidates Who Question Elections",
            Variable == "norm_loyalty_perception" ~ "Perceived Out-Party Belief that Out-Partisans Should\nbe Loyal to Out-Party Candidates Who Question Elections",
            Variable == "norm_polling" ~ "Belief that In-Party Should\nReduce Polling Stations in Out-Party Areas",
            Variable == "norm_polling_perception" ~ "Perceived Out-Party Belief that Out-Party\nShould Reduce Polling Stations in In-Party Areas",
            Variable == "out_party_therm" ~ "Warmth\nToward Out-Party",
            Variable == "violence_protest" ~ "Support for In-Partisan Who Led\nan Illegal Protest Against the Out-Party",
            Variable == "violence_vandalism" ~ "Support for In-Partisan\nWho Vandalized Out-Party Signs",
            Variable == "violence_assault" ~ "Support for In-Partisan Who\nThrew Rocks at Out-Partisans",
            Variable == "violence_assault_perception" ~ "Perceived Out-Party Support for\nOut-Partisan Who Threw Rocks at In-Partisans",
            Variable == "violence_arson" ~ "Support for In-Partisan Who\nSet Fire to Out-Party Headquarters",
            Variable == "violence_assault_deadly_weapon" ~ "Support for In-Partisan\nWho Drove into Out-Party Crowd",
            Variable == "violence_murder" ~ "Support for In-Partisan Who\nStabbed Prominent Out-Partisan",
            Variable == "violence_murder_perception" ~ "Perceived Out-Party Support for\nIn-Partisan Who Stabbed Prominent Out-Partisan"
        )
    ) %>%
    arrange(desc(estimate))

# Factorize variable column
boot_results$Variable <- boot_results$Variable %>%
    factor(
        levels = boot_results$Variable
    )

# Identify highlighted sections
boot_results <- boot_results %>%
    mutate(
        rect_left = case_when(
            grepl("^Perceived", Variable) ~ min(lower)
        ),
        rect_right = case_when(
            grepl("^Perceived", Variable) ~ max(upper)
        ),
        rect_bottom = case_when(
            grepl("^Perceived", Variable) ~ as.numeric(Variable) - 0.5
        ),
        rect_top = case_when(
            grepl("^Perceived", Variable) ~ as.numeric(Variable) + 0.5
        )
    )

# Calculate summary statistics
boot_results %>%
    subset(grepl("^Perceived", Variable)) %>%
    pull(estimate) %>%
    range() %>%
    round(2)

boot_results %>%
    subset(!grepl("^Perceived", Variable)) %>%
    pull(estimate) %>%
    range() %>%
    round(2)

############
### PLOT ###
############

boot_results %>%
    ggplot() +
    geom_pointrange( # Points have to be plot first for rectangles to work
        aes(
            x = estimate,
            xmin = lower,
            xmax = upper,
            y = Variable,
            color = grepl("^Perceived", Variable)
        )
    ) +
    geom_rect(
        aes(
            xmin = rect_left,
            xmax = 26,
            ymin = as.numeric(Variable) - 0.5,
            ymax = as.numeric(Variable) + 0.5
        ),
        fill = "grey",
        alpha = 0.2
    ) +
    geom_pointrange( # Duplicate points so not obscured by rectangles
        aes(
            x = estimate,
            xmin = lower,
            xmax = upper,
            y = Variable,
            color = grepl("^Perceived", Variable)
        )
    ) +
    geom_text(
        data = global_means,
        aes(
            y = Variable,
            x = max(boot_results$rect_right, na.rm = TRUE) + 1.5,
            label = round(numeric.mean, 1)
        ), size = 3
    ) +
    annotate(
        "text",
        x = max(boot_results$rect_right, na.rm = TRUE) + 1.5,
        y = 22.5,
        label = "Global Mean",
        size = 3,
        fontface = "bold"
    ) +
    scale_y_discrete(expand = expansion(add = c(0, 1))) +
    scale_color_manual(values = c("black", "#56B4E9")) +
    theme_prl() +
    theme(
        axis.title.x = element_markdown(),
        axis.title.y = element_markdown(),
        axis.text.y = element_text(size = 8)
    ) +
    labs(
        x = "**Over-Time Variability in Individuals' Responses**<br>(Average Within-Individual Standard Deviation as % of Outcome's Range)",
        y = "**Variable**"
    ) +
    guides(
        color = "none",
        size = "none"
    )

# Save plot
ggsave(
    "Images/within_ind_SDs.png",
    width = 12.3, height = 9
)

#############
### TABLE ###
#############

boot_results %>%
    select(Variable, estimate, lower, upper) %>%
    arrange(estimate) %>%
    mutate(
        Variable = Variable %>% str_replace("\\\n", " ")
    ) %>%
    mutate(
        `Average Within-Individual SD` = paste0(
            estimate %>% round(2) %>% as.character(),
            " (",
            lower %>% round(2) %>% as.character(),
            ", ",
            upper %>% round(2) %>% as.character(),
            ")"
        )
    ) %>%
    select(Variable, `Average Within-Individual SD`) %>%
    datasummary_df(
        title = "Estimates from Figure 1",
        notes = "Note: Average within-individual standard deviation for each
        variable. Parentheses indicate the bootstrapped 95% CI of each estimate
        (based on 2,000 simulations). For each variable, we calculated the
        standard deviation of each panelist's responses (across panel waves)
        and then averaged these individual-level standard deviations.
        All variables were re-scaled (0--100) to ensure comparability. To
        account for varying durations between different panelists' responses,
        in calculating averages, individuals' standard deviations were weighted
        by $1/x$, where $x$ is the standard deviation of the individual's
        survey response dates.",
        output = "Tables/within_individual_SDs.txt"
    )
